home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 105_01 / ctrig.c < prev    next >
Text File  |  1984-06-03  |  6KB  |  242 lines
  1.  
  2.  
  3.  /*  **** CTRIG.C *****
  4. A group of programs in C, using the BDS-C floating point package,
  5. as modified by L. C. Calhoun called FLOATXT, which compute some
  6. commonly used transcendental functions - to wit*
  7.     sine, cosine, tangent and arctangent
  8.     convert degrees to radians, convert radians to degrees
  9. These functions are discussed in detail in CTRIG.DOC
  10.  
  11. L. C. Calhoun
  12. 257 South Broadway
  13. Lebanon, Ohio 45036   513 932 4541/433 7510
  14.  
  15.   */
  16.  
  17. /* simple ones first converting degrees - radians */
  18.  
  19. char *degtorad(rad,deg) /*obvious arguments in 5 char fp */
  20. char *rad, *deg;
  21. {
  22.     char *fpmult(),radindeg[5];
  23.     initb (radindeg,"217,27,125,71,251");
  24.     fpmult(rad,deg,radindeg);
  25.     return (rad);
  26. }
  27.  
  28. char *radtodeg(deg,rad) /* 5 char fp arguments */
  29. char *deg, *rad;
  30. {
  31.     char *fpmult(),deginrad[5];
  32.     initb (deginrad,"10,114,151,114,6");
  33.     fpmult(deg,rad,deginrad);
  34.     return (deg);
  35. }
  36.  
  37. /* service function sinev which evaluates when range of angle
  38. reduced to plus or minus pi/2 (90 deg) */
  39.  
  40. char *sinev(result,angle)
  41.  
  42. char *result, angle[5];
  43. {
  44.     char *fpmult(),x[5],xsq[5];
  45.     char coef[5][5],termreslt[5];
  46.     char *fpadd(),*fpasg();
  47.     int index;
  48.  
  49.     /*  use the exponent part of the floating point
  50.         to check for threat of underflow  use small
  51.         angle approximation if appropriate  */
  52.     if ( (angle[4] > 128) && (angle[4] < 226) )
  53.        {fpasg(result,angle);       return (result);
  54.       }   /* solution to fpmult underflow problem */
  55.  
  56. /* series coef are 1., -.1666666, .008333026, -.0001980742,
  57.     .000002601887  determined from coefset program */
  58.     initb(coef[0],"0,0,0,64,1");
  59.     initb(coef[1],"111,170,170,170,254");
  60.     initb(coef[2],"185,162,67,68,250");
  61.     initb(coef[3],"208,234,38,152,244");
  62.     initb(coef[4],"166,13,78,87,238");
  63.     fpasg(x,angle);
  64.     fpmult(xsq,x,x);
  65.     setmem(result,5,0);
  66.  /* to this point the coef have been initialized, the angle
  67.     copied to x, x squared computed, and the result initialized */
  68.  
  69. /* now to do the polynomial approximation */
  70.     index = 0;
  71.     while ( (index <= 4) && ( (x[4] > 194) || (x[4] < 64) ) )
  72.     /* use index for loop, and exponent of x to avoid underflow
  73.        problems */
  74.     {fpmult(termreslt,coef[index],x);
  75.     fpadd(result,result,termreslt);
  76.      index++;
  77.      fpmult(x,x,xsq);
  78.     }
  79.     return (result);
  80. }
  81.  
  82.  /* here is sine(result,angle) with angle in radians */
  83.  
  84. char *sine(result,angle)
  85.  
  86. char *result, *angle;
  87. {
  88.     char *fpmult(),twopi[5],halfpi[5];
  89.     char mtwopi[5],mhalfpi[5],*fpasg(),*fpchs();
  90.     char pi[5],*sinev(),*fpadd();
  91.     char y[5],*fpsub();
  92.     int fpcomp(), compar;
  93.     int signsine;
  94.  /* some initialization required here */
  95.     signsine = 1;
  96.     initb(twopi,"121,238,135,100,3");
  97.     initb(halfpi,"121,238,135,100,1");
  98.     initb(pi,"244,239,135,100,2");
  99.     initb(mtwopi,"135,17,120,155,3");
  100.     initb(mhalfpi,"135,17,120,155,1");
  101.     fpasg(y,angle);
  102.     while (fpcomp(y,twopi) >= 0)
  103.        {fpsub(y,y,twopi);
  104.        }
  105.     while (fpcomp(y,mtwopi)<= 0)
  106.        {fpadd(y,y,twopi);
  107.        }
  108.     if(fpcomp(y,halfpi) > 0)
  109.        {signsine *=-1; fpsub (y,y,pi);
  110.        }
  111.     if(fpcomp(y,mhalfpi)<  0)
  112.        {signsine *=-1; fpadd (y,y,pi);
  113.        }
  114.     sinev(result,y);
  115.     if (signsine > 0) return (result);
  116.  /* minus so need to change sign */
  117.     else return ( fpchs(result,result) );
  118. }
  119.  
  120.  /* cosine(result,angle) with angle in radians  - uses sine */
  121.  
  122. char *cosine(result,angle)
  123.  
  124. char *result, *angle;
  125. {
  126.     char *sine(),*fpsub(),halfpi[5],y[5];
  127.     initb(halfpi,"121,238,135,100,1");
  128.     fpsub(y,halfpi,angle);
  129.     sine(result,y);
  130.     return (result);
  131. }
  132.  
  133. /* tangent(result,angle) with angle in radians, returns very 
  134. large number for divide by zero condition */
  135.  
  136. char *tangent(result,angle)
  137. char *result, angle[5];
  138. {
  139.     char *sine(), *cosine(), *fpdiv(), zero[5];
  140.     char sresult[5], cresult[5], intres[5], big[5];
  141.     char *fpmult(), *fpmag();
  142.  
  143.     sine(sresult,angle);
  144.     cosine(cresult,angle);
  145.     /* check magnitude of denominator :*/
  146.     /* check magnitude of denominator using exponent */
  147.     if ( (cresult[4] > 128) && (cresult[4] < 132) )
  148.        {initb(big,"30,207,228,127,128"); /*big number */
  149.         if ( sresult[3] > 127 ) /*use mantissa sign */
  150.            return ( fpchs(result,big) );
  151.         else return ( fpasg(result,big) );
  152.        }
  153.     /* check for small angle, use small angle approx to
  154.        avoid underflow */
  155.     if ( (angle[4] < 226) && (angle[4] > 128) )
  156.        return ( fpasg(result,angle) );
  157.     else
  158.        return ( fpdiv(result,sresult,cresult) );
  159. }
  160.  
  161. /* atanev(result,x) evaluates arctangent for 0 <= x < infinity,
  162.     result in radians */
  163.  
  164. char *atanev(result,x)
  165.  
  166. char  result[5], x[5];
  167. {
  168.     char coef[8][5],piover4[5],y[5];
  169.     int index;
  170.     char *fpadd(),*fpmult(),*atof(),ysq[5],one[5];
  171.     char yterm[5],termreslt[5],*fpasg();
  172.  
  173.     initb(piover4,"121,238,135,100,0");
  174.     initb(one,"0,0,0,64,1");
  175.  
  176.     /* small angle approximation */
  177.     if ( (x[4] > 128) && (x[4] < 226) )
  178.     /* use fp exponent to check size, use small angle */
  179.      return ( fpasg(result,x) );
  180.     else
  181.     fpasg(result,piover4);
  182.  
  183.     /* check for argument near one */
  184.     fpsub(yterm,x,one);
  185.     /* if it is close to one, as checked by exponent,
  186.        return the result pi/4  */
  187.     if ( (yterm[4] > 128) && (yterm[4] < 243) )
  188.        return (result);
  189.  
  190.  
  191.     initb(coef[0],"184,254,255,127,0");
  192.     initb(coef[1],"191,239,172,170,255");
  193.     initb(coef[2],"70,86,32,102,254");
  194.     initb(coef[3],"102,203,201,184,254");
  195.     initb(coef[4],"28,240,187,98,253");
  196.     initb(coef[5],"134,24,127,141,252");
  197.     initb(coef[6],"108,44,139,89,251");
  198.     initb(coef[7],"55,10,148,189,249");
  199.     fpadd(termreslt,x,one);
  200.     fpdiv(y,yterm,termreslt);
  201.  
  202.     fpmult(ysq,y,y);
  203. /* do poly evaluation */
  204.     index = 0;
  205.     /* poly evaluation checked by index limit, and check of
  206.     variables for under/over flow */
  207.     while ( (index <= 7) && ( (y[4]<100) || (y[4]>140) ) )
  208.      {fpmult(termreslt,coef[index],y);
  209.       fpadd(result,result,termreslt);
  210.       index++;
  211.       fpmult(y,y,ysq);
  212.      }
  213.  
  214.     return (result);
  215. }
  216.  
  217. /* arctan(result,angle) is floating point arctangent evaluation */
  218.  
  219. arctan(result,x)
  220. char result[5], x[5];
  221.  
  222. {
  223.     char *atanev(),*fpasg(),y[5];
  224.     char *fpchs(),halfpi[5];
  225.     int index;
  226.     /* check exponent for very large argument */
  227.     if ( (x[4] > 100) && (x[4] <= 128) )
  228.        {initb(halfpi,"121,238,135,100,1");
  229.         fpasg(result,halfpi);
  230.        }
  231.     else  /* go through evaluation */
  232.       {fpmag(y,x);
  233.        atanev(result,y);
  234.       }
  235.  
  236.     if ( x[3] > 127 )
  237.       {return ( fpchs(result,result) );}
  238.     else return (result);
  239. }
  240.  
  241.  
  242.